We study a bilevel \emph{max-max} optimization framework for principal-agent contract design, in which a principal chooses incentives to maximize utility while anticipating the agent's best response. This problem, central to moral hazard and contract theory, underlies applications ranging from market design to delegated portfolio management, hedge fund fee structures, and executive compensation. While linear-quadratic models such as Holmstr"om-Milgrom admit closed-form solutions, realistic environments with nonlinear utilities, stochastic dynamics, or high-dimensional actions generally do not. We introduce a generic algorithmic framework that removes this reliance on closed forms. Our method adapts modern machine learning techniques for bilevel optimization -- using implicit differentiation with conjugate gradients (CG) -- to compute hypergradients efficiently through Hessian-vector products, without ever forming or inverting Hessians. In benchmark CARA-Normal (Constant Absolute Risk Aversion with Gaussian distribution of uncertainty) environments, the approach recovers known analytical optima and converges reliably from random initialization. More broadly, because it is matrix-free, variance-reduced, and problem-agnostic, the framework extends naturally to complex nonlinear contracts where closed-form solutions are unavailable, such as sigmoidal wage schedules (logistic pay), relative-performance/tournament compensation with common shocks, multi-task contracts with vector actions and heterogeneous noise, and CARA-Poisson count models with $\mathbb{E}[X\mid a]=e^{a}$. This provides a new computational tool for contract design, enabling systematic study of models that have remained analytically intractable.

In this paper, we present a novel general framework grounded in the factor graph theory to solve kinematic and dynamic problems for multi-body systems. Although the motion of multi-body systems is considered to be a well-studied problem and various methods have been proposed for its solution, a unified approach providing an intuitive interpretation is still pursued. We describe how to build factor graphs to model and simulate multibody systems using both, independent and dependent coordinates. Then, batch optimization or a fixed-lag-smoother can be applied to solve the underlying optimization problem that results in a highly-sparse nonlinear minimization problem. The proposed framework has been tested in extensive simulations and validated against a commercial multibody software. We release a reference implementation as an open-source C++ library, based on the GTSAM framework, a well-known estimation library. Simulations of forward and inverse dynamics are presented, showing comparable accuracy with classical approaches. The proposed factor graph-based framework has the potential to be integrated into applications related with motion estimation and parameter identification of complex mechanical systems, ranging from mechanisms to vehicles, or robot manipulators.